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Abstract 



An explicit numerical scheme is proposed for solving the initial-boundary value problem for 
the radiative transport equation in a rectangular domain with completely absorbing boundary 
condition. An upwind finite difference approximation is applied to the differential terms of 
the equation, and the composite trapezoidal rule to the scattering integral. The main results 
' are positivity, stability, and convergence of the scheme. It is also shown that the scheme can 

be regarded as an iterative method for finding numerical solutions to the stationary transport 
equation. Some numerical examples for the two-dimensional problems are given. 



?3 ■ 1 Introduction 



The radiative transport equation (RTE) is an integro-differential equation which is widely used 
as a mathematical model of the flow of particles traveling through a medium which absorbs and 
scatters the particles [3] . For example, let us regard the propagation of light in a turbid medium 
as a flow of photons and let I = I(t, x, £) denote the density of the photons moving in the direction 
£ at time t and position x. With suitable assumptions on the optical property of the medium, we 
can obtain an equation of the form (|la[) (given below) which governs the time evolution of I. 

Let fl be a rectangular domain with edges parallel to the coordinate axes in the d-dimensional 
Euclidean space M. d (d — 2 or 3) and let S^ 1 denote the unit sphere in M. d . We consider the 
£f*^ . initial-boundary value problem 

-7-7757 = -£ • V X J - {fi a (x) + fi s (x)}I + n s (x) f p(x; £, £')/(*, x, £') du e + q(t, x, £), 

C{X) Ol J gd-1 

X ■ < t < T, (x, e Q x S d -\ (la) 

1(0, x, o = i (x, 0, (x, € n x s d -\ (ib) 

l(t,x,0=h(t,x,0, 0<(<T,(i,fleL, (lc) 

where 

T_ := {(£,£) efflx S 4 - 1 ; n(x)-(< 0} , 

with n(x) the outward unit normal vector at x € dQ. In equation (jlap the dot ■ indicates the 
standard inner product in R d , V x denotes the spatial gradient and doj^> the surface element on 
The coefficient c(x) denotes the speed of particles at position x, fia,(x) the absorption 
coefficient, and fJ, s (x) the scattering coefficient. The integral kernel p(x\ £, £') is called the scattering 
phase function and expresses the density of the conditional probability of the event that a particle 
moving in the direction £' is scattered into the direction £ at position x. The inhomogeneous term 
q(t,x,£) is the density of internal source of particles emitted in the direction £ € at time t 

and position x. 

An application of the RTE of the form (fTa|) can be found in researches on diffuse optical 
tomography (DOT), which is expected to be a new modality of noninvasive medical imaging which 
uses low-energy light [TJ [2] . The propagation of light in the biological tissue under examination is 
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modeled by equation ([lap whose coefficients are unknown quantities to be determined by DOT. 
The coefficients represent optical properties of the biological tissue, and they actually depend not 
only on position x but also on time t. However, the speed of light is so faster than time scale 
of tissue activity that we neglect their dependence on t. The procedure of imaging by DOT is 
regarded as an inverse problem of recovering the unknown coefficients by observing the solution / 
at the boundary of the domain. When we analyze the inverse problem, it is extremely valuable to 
have an efficient numerical method for solving the problem (JTJ since it is practically impossible to 
have a simple analytic expression for the solution of ([1]). 

The purpose of the present paper is to construct an explicit numerical scheme for solving 
the initial-boundary value problem (JTJ) - In our discretization of (|la[) . the differential operator 
c~ 1 dt + £ ■ V 2 is approximated by an upwind finite difference, and the scattering integral by the 
composite trapezoidal rule. We describe the detail only for the two-dimensional case in Sj2j and give 
a brief account of the three-dimensional case in S}5J The main results are the stability (Theorem 
2.21 Theorem 15. ip and the convergence (Theorem 12.41 Theorem 15. 2p of the scheme. In ^3]we also 
show that the scheme can be regarded as an iterative method for finding numerical solutions to the 
stationary transport equation. In <T4]we give examples of numerical solution of the initial-boundary 
value problem ([1]) and the corresponding stationary problem in the two-dimensional case. 

Finally we formulate our assumptions on the data of the problem ([1]) which are fixed throughout 
the paper. We assume that the coefficients c, /i a , and fi s are measurable functions satisfying 



< c(x), 0<H a (x), 0<Li s (x), n s (x)^0 1 
+ — c^^A ^ „* — i n f 



c + := sup c(x) < +oo 



> 0. 



{x- M*)#0} Hg(x) 

By its physical interpretation, the scattering kernel p(x;£, £') must be nonnegative and satisfy 



(2) 



We assume moreover that p depends on £ and £' through the angle, say a, formed by them, so 
that we can write 

2?(s;£)0 =p(x;a), (3) 

with a nonnegative function p(x; a) which is 27r-periodic and even in a for each x\ we also call 
p(x; a) the phase function. Consequently, the condition <j2j) is equivalent to 
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and is further rewritten in terms of p as 

,-2tt 



/■ Z7T p 7T 

/ p~(x;a)da = l for d = 2; 2ir / p(x; a) sin a da = 1 for d — 3. 
Jo Jo 



We assume that p(x; a) is of class C 2 in a for each x € fl and that 

d 2 p 



da 2 



:= sup 

oo (x,a)eOx [0,2ir] 



d 2 p 



< oo. 



The inputs are assumed to be real bounded measurable functions: 

q(t,x,0 e L°°((0,T) x n x s 11 - 1 ), l (x,0 e x S^ 1 ), h(t,x,0 e L°°((0,T) x r_). 

Their L°°-norms are denoted by ||q|joo, ll^olloo, and ||/i||oo- 
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2 Upwind finite difference scheme for RTE and its proper- 
ties for the two-dimensional case 

Let fl = (0,Li) x (0,L 2 ) be a rectangle in R 2 . In this section and the next we put X — Q x S 1 , 
which is the phase space of equation (fTaj) . For three positive integers Mi, M 2 , and M and a 
positive number At, we put 

Ax 1 =L 1 /M 1 , Ax 2 =L 2 /M 2 , A8 = 2tt/M, 

and consider the grid points (tfc,a:«,£n) defined by 

t fe = fcAt, ■ = (iAxi,jAx 2 ), £ n = (£n,i> Cn,2) = (cosnA0, sinnA0), 

where k,i,j,n G Z. Using as the value corresponding to I(tk,Xij,£ n ), we discretize the 

problem ([TJ as follows: 



, jk+l _ jk 

c(xij ) At 



,,7, ft ~^ Qi,j } n ' 





J*. 



0<fc<T/At-l, (i^Qel, 
■ h(t k ,Xij,£ n ), < ft < T/At, (ly^JeL, 



where n = g(t/t, and the operators ^4a, Ea, and if a are defined by 



4 Tfc _ t I i+l,j,n I i-l,j,n , 1* 
^W,\7, „ - r- K 



2Axi 

jk jk 



Sn,l 



rfc o r A; 1 jk 

i-\-l,j,n i.j,n * *i—l,j,n 

2Axi 



; — 2I k ->- I 
£n,2 — " + Cn,2 — ' 



2Ax 2 

Ea^j> = {/i s (xy) + /i a (Xy-)}7j j 



2Ax9 



M-l 



(4a) 
(4b) 
(4c) 



(5a) 
(5b) 

(5c) 



Remark 2.1. (i) The operator Aa is defined so that the combination 



J- i,j,n i : j,n 

c(xij) At 



(6) 



gives the first-order upwind finite difference approximation to the partial differential operator 
c~~ 1 d t + £ • V x - For example, for n with £ n> x > and £ n , 2 > 0, ([6]) is equal to 



c(xij ) At 



i.j,n ^i.j^n j. ^i.j^n i—lj.n j. i.j,n i,j—X,n 
Sn,l — ■ a r?„,2- 



Axi 



Ax 2 



The other three cases can be checked similarly, 
(ii) Scheme (j4]) is explicit; see (fTT]) . 

Theorem 2.2. (Stability and positivity for d — 2) Suppose that the discretization parameters 
satisfy 



c+At c+At 

■ + -: < 1, 



Axi Ax 2 



d 2 p 
da 2 



A6 2 <- 



(7a) 
(7b) 
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Then we have 

ll^lloo < ||/o|| oo H~ 1 1 ^1 1 1 oo H~ c + ||g||ooT, 0<fc<T/Ai, 

where 

Halloo := sup {\I^ n \; fegauL}. 
Moreover if q > 0, I > 0, and I± > everywhere, then ijv n > holds for all fc, i, j, n. 

In the proofs below we write c, /^ a , and fi s instead of c(xij), /j, a (xij), and fi s (xij), respectively. 
Proof. We show that 

Halloo < Poll c + \\q\\oot k (8) 

for < k < T/At by induction on k. It is obvious for k — 0. Assume now that (|8|) holds for some 
k < T/At — 1. The proof shall be finished by showing that 

|^+n| ^ ll^lloc + c + ||g||ooAi, (zii,£ n ) G A:. (9) 
Indeed, it follows from (jHJ and (|9|) that 

< \\IoWoo + ||/i||oo + c + ||g|| 0O t fc +i, (lij^Jel 
Since = |/i(f/e, ly, £ n )| < Pi||oo for 0%'>£n) G r_, we obtain 

||/ fc+1 ||oo < INloo + ||/l||oo + C+Hgllooi/fe+l. 

We now show ©. We introduce an auxiliary operator by the relation 



or explicitly 



n rfe (ie,.ii-^,i)^+i.,,„ + (l6Ml+c«,i)^-i 



2Axi 

(|£ra,2| - Cn,2)-fj,j+l, n + (|£n,2| + ^n.a)^-!, 



2Ax 2 

Then we write (|4aj) as 

rfe+l 



{l + cAi( Ms+Ma )}/^ 
I— -I 



We find from (Taj that 

, * , cAi . „ . cAt 

i-IWI^-l^ 2 |^>o, 

and from (jTUJ) that 

rfe 



7> ,/ | I lC",l l . |€tX,2 I I Mr/, 

Ba/«.„|< l^ + ^lll-f 



I oo 1 



hence from (jTTJ) that 



{l + cAi^+Ma)}!/^, 



(10) 



1 " |£»,1 I XT - IWIXT ) + cAi ( BAl Ln + K ^Ln + Qi,j,n) ■ (11) 



< Wl^+cAtlK^J+cAtWqlU (12) 
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By (|5c|) and ([3]) we have 

M-l M-l 

I^aJ^I < m,A0 ^ p(x«;e»,^)||/ fc ||oo = H^Hoo^Afl P( x ij;0» ~ *«)• (13) 

t/=0 i/=0 

Here we use the next lemma. 

Lemma 2.3. Let f(9) be a 27r-periodic function of class C 2 . Let M be a positive integer and 

A0 = ^, 9 m = mA9, m = 0,1,..., Af-l. 
Then there exists a point 6*' € [0,27r) such that 

Af— l 



/ * d9 - A9J2 Wm) = ^/"(^')A0 2 

^° m=0 



This can be shown routinely by taking the periodicity of / into account. 
Applying Lemma 12.31 to the function f(9) — p(xij\ 9 — 9 n ), we find 



p2TT M-l q 2 

/ p{x if ,e - 9 n ) d9~A9Y] p(ssy; 0„ - n ) = — A9 2 ^p{x^9 - 9 n ) 
Jo v=0 U CIO 

for some 9' = 9'ij, n - Since the integral of p is equal to 1, we get 



M-l ~ 2 

A# £ pfc,-; 0„ - n ) = 1 - ^A0 2 Q^pixijie - 0„) 

u=0 



<1 + M* (14) 



by using (I7bp . Combining (fl^j) — (p~¥)) and observing that /x s /x* < p, a on f2, we obtain 
{1 + cAi(/i s + n*)} \ < {1 + cAt(^ s + ^a)}||/ fe ||oo + cllffHooAt. 

Dividing both sides by 1 + cAt(/i s + jtx a ), we arrive at (J9)). 

The positivity is shown as follows: if lfj n > for all i,j, n at a time level k, then we have 
/\ a/;;, „ > by (5cJ), and B A I^ jn > by ^TUJ) . Hence > at the next time level fc + 1 by 

CP and g > 0. " " □ 

Theorem 2.4. (Convergence for d = 2) In addition to the hypotheses of Theorem 12.21 we 
suppose the following three conditions: (i) the mesh ratios At/Axi and At/Ax2 are kept fixed; 
(ii) there exist positive numbers p,J and fij such that 

< /i a (x) < jti+, < /i s (:r) < a; G 0; 

and (hi) the solution I(t,x,£) of |T]) belongs to C 2 ([0,T) x (X U L_)) and is bounded together 
with its derivatives of the first and the second order. Then we have 

||/(*fc, •, ') - I k \\oo < C(At + A9 2 ), < k < T/ At, 
where C is a positive number independent of At and A9. 

Proof. Let I(t, x, £) be the solution of (p} satisfying (iii). We define the local truncation error rf j n 
by the relation 

1 lipk+l i Xij , £ n ) -f , Xjj , £n) 

c At 

= AA/(*fe,Xij,$ n ) - ^Al(tk+l,Xij,£ n ) + K A I{tk,Xij,^ n ) + 9ij> + T i,j,m (I 5 ) 
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where the operators Aa, Saj and Ka are applied to I(tk,Xij,£ n ) in place of By using 

equation (fTaj) . Taylor's theorem, and by applying Lemma 12751 to the function 

f{6) = p(xij;£ n ,r))I(t k ,Xij,ri) with rj = (cos 6, sin 0), 

we arrive at 



19/,, x . £ n \ d I , . j. \ a tin 2 & I f n j~ \ a 



2 9xi 2 
91 - ^A/ia 



2 <9x 2 2 



(Ms + ^)-^(tk,^3^n)At + -jyA6> —p(x ij ;^n,v)I(tk,x ij ,r]) 

±1 +» n „A „ll 



for some i' fc , i' fc ', xL, and x"-. Multiplying this equality by c and using the assumptions, we obtain 



c+ sup |r* ? J<C"(At + A0 2 



k,i,j,n 



(16) 



with C" independent of At and A#. 

Next we consider the discretization error J5* ■ n := 7^ - n — I(tk,Xij,£ n ). We subtract (fT5j) from 
(l4a| to have 

~ ^ — ^i^tin — A ij,n + A A fl i,j,n — T i,i,"' 

o<fc<r/Af-i, 

We observe that £7 = {-E^. n } satisfies equation (Ha]) with gfj- n = — Tijm and the homogeneous 
initial-boundary conditions. Therefore the proof of Theorem 12.21 shows that 



\E 



k+l I 



Then P^|) gives 



with C := TC 



< + + Tc+ sup \T? j n \, < k < T/At - 1. 



l^lloo < C(At + A6» 2 ), 0<fc<T/At, 



□ 



Remark 2.5. Here we consider a particular situation: p(x; 0) is independent of x and analytic in 
with period 27r. We write p{9) = p(x; 0) and expand it into the Fourier series 



p(9)= J2 c ^ me 



n— — 00 



2^ Jo 



p(d)e~ ln9 d9 



(17) 



Since p{9) is analytic, there exist numbers C > and < r < 1 such that (see [8]) 

|c„|<CY |n| , neZ, 



/ ji(0)d0 - P(6m) 

J u=0 



< 



AnCr 



M 



1 - r 



M 



M G N, A9 



2tt 



Then the conclusions of Theorem 12.21 and Theorem 12.41 can be obtained if we assume 

4ttCV m 

< Mi 



1 - r 



M 



(18) 



(19) 



instead of (ffb|). Indeed, since / p(0) dO = 1, we find from (JISJ and (JTSJ) that 

Jo 

M-l 

A0 ^ p(0„ - 



<1 + M*, 



u=0 



which corresponds to (|14j) . The proofs of the theorems are valid without further modification. 
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Example 2.6. In researches on diffuse optical tomography (TJ [2J, the main interest is in the 
three-dimensional case, and the Poisson kernel 

P( g ) = J- n o 2W2 ' with °<5<1, 

47T (1 — 2g cos # + g z yi z 

is frequently used as the phase function. The Poisson kernel is often referred to as the Henyey- 
Greenstein kernel in that field and the parameter g indicates the strength of forward scattering by 
the medium [5]. If g = 0, then p(9) — 1/Air identically and the particles will be scattered uniformly 
in all directions. If g gets closer to 1, the proportion of photons scattered in directions near the 
incident one gets larger. 

Since we are dealing with the two-dimensional case in this section, we consider the two- 
dimensional analogue: 



1 1 



■2 



6(6) = ~ with < g < 1. 

Fy ' 2?r 1 - 2gcos<9 + g 2 ~ J 

Its Fourier series expansion is 

PW = 5- E 9 ]n ^ n6 , i-e., c„ = ^— in (HZ}, 

Z7T ' — ' Z7T 

n— — oo 

so that, by Remark I2.51 the scheme is stable if Axi, Ax2, At, and M satisfy (jTaj) and 

2o M 

Y^m < M*. (20) 
If <7 = 0, in particular, the condition (j2TJ)l trivially holds, so that the stability is independent of M. 



3 Numerical solution of the stationary problem for RTE 



We turn our attention to the initial-boundary value problem ([TJ for an infinite interval of time 
< t < oo with time-independent inputs q — q(x,£) and 1% = Ii(x, £). We consider the scheme 
(HI for all k = 0, 1, 2, . . . making an additional assumption that 

(c^ a )~ := inf c(x)fj, a (x) > 0. 

Our aim here is to show that the scheme ^ serves an iterative method for approximating the 
solution of the corresponding stationary problem: 

V X J— (fi s (x) + fi a (x)) J + » s (x) [ p(x; C, ?)J(x, O fa? + = 0, (x, £) e X, (21a) 

Js 1 



(21b) 

(22a) 
(22b) 

with Ji,j >n corresponding to J(xij,£ n ). Then (|22|) is uniquely solvable for sufficiently small A(9 as 
the next theorem shows. 



j(x,0 = h(x,0, (i,0er_, 

We discretize the problem (|2"Tj) as follows: 

(A A - E A + K A )J itjin + q iij>n = 0, (x.y , £ n ) € X, 



Theorem 3.1. Suppose that A6 satisfies 

||p||ooA0 < 1 and 



d 2 p 



d9 2 



A9 2 < —p* 

7T 



Then there exists a unique solution Ja = {J%,j,n} °f (|22p and it satisfies 



|| Ja||oo := SUp \Ji,j, n \ < ll^llloo + C||g| 
i,j,n 

where C is a positive number independent of Ax±, Ax2, and A9. 



(23) 



7 



We refer to [I] for a proof of Theorem 13.11 and instead show that the scheme (0| generates an 
iterative method for solving 



Theorem 3.2. Suppose that the discretization parameters for the scheme ((4]) are fixed so that 



c+At c+At 



< 1, ||p||ooA0 < 1, and 



dp 



da 2 



Ax i Ax 2 

Then there exists a number p, < p < 1, such that 

\\I k - Ja||oo < (||/o||oo + ll'l||oo)/, k > 



AO 2 < —fi* 

7T 



Proof. We put R$ in = A n - J iJtn . Writing (gSaJ) as 



~ c At 



(A A -S A +ifA)Ji,. 



and subtracting it from (I4al) , we have 



fe+i 



Let A be a number, < A < 1, determined by 



d 2 p 



da 2 



A9 2 = 

7T 



Reasoning as in the proof of Theorem 12.21 we obtain 

>fc+i 



{l + cAt(fi s + ^)}\Rf+ L n \ <{l + cAt(rt, + A/x»)}||i2*||oo, fc>0, (i l3 ,C„)eX 
Since RH^ = for (;£„,£„) £ T_, we have 



\R 



<( S u P ;+^;+ Ac ^f i ii/?^ 



Observing that /i s < /i a //i* and c/i a > (c/x a ) on J7, we see that 

1 + c^ s Ai + Xc^At < 1 + cfi a At/n* + \cfi a At _ 1 + cfi a At(l/fi* + A) 
1 + c/i s At + c/i a Ai ~ 1 + c^At/^i* + c^ a Ai ~ 1 + cfx a At(l/fi* + 1) 
l + (c^ a )-At(l//i* + A) 



< 



p<l. 



1 + (cA* a )-A*(l//x* + 1) 
Thus we obtain ||i? fe+1 ||oo < Hl^lloo for all fc > and hence 

||/ fc - JaUoo = ||i? fc ||oo < /Halloo = P k \\l° - JaUoo, k > 0. 

The proof completes by using (|23|) . 



□ 



4 Numerical Examples 

Here we consider a square domain fl = (0, 50) x (0, 50) C M 2 of sides 50 mm and assume that the 
domain is occupied by a medium of uniform optical properties: /i a = 0.08 mm , /i s = 1.09 mm , 
and c = 0.196 mm/ps. Units of length and time in our numerical computation are chosen to be 
millimeter and picosecond, respectively. We employ the Poisson kernel 



m 



i 



2-7T 1 — 2g cos 9 + g 



3 = 0.9 
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as the phase function. The boundary condition is set to 



' exp (- ( 6 I , 24.9 <x 1 < 25.1, x 2 = 0, 



0, otherwise, 

where a — 0.2 and £ = (cos 9, sin 9). The initial condition is Iq(x,^) = in X, and Io(x,£,) and 
I\{t, x, £) satisfy compatibility on T_. The parameters and the inputs are chosen so as to simulate 
a phantom experiment carried out at Human Brain Research Center, Kyoto University [7]. 

Numerical computation is processed with At = 0.1, Ax = Ay = 0.1, and A9 = 27r/60. These 
satisfy the stability conditions (|7aj) and (|20|) . since 



60 = A/ > log ^ ~ 31.71. 

Fig. [1] shows the numerical solution /* 3 - n for several points iy. The solution is shown in the polar 
coordinate with respect to 9 for £ = (cos 0, sin 0). Fig. [2] shows contour lines of integrated light 

intensity / I(t,x,£)dw£ at time t and position x. 
Js 1 

As we have seen in fj3[ the sequence of grid functions {/ fc }^ generated by the scheme ([4]) 
converges to the solution Ja of (|22|) as /c tends to infinity. The convergence of {/ fc } to Ja is shown 
in Fig. [3l in which the horizontal axis is the number of time step k and the vertical axis is the 
difference in logarithmic scale. The cross signs (x) show ||/ fe — I fc_1 ||oo and the plus signs (+) show 
||/ fe — Ja||oo- Both are approximately proportional to 0.99807 fc , which means that {I k } converges 
to Ja exponentially fast. The difference — Ja||oo is saturated around 10 -12 since the solution 
Ja of (f2"2"j) is found by an iterative method with tolerance 10~ 12 . On the other hand, \\I k — / 1 ||o Q 
is not less than 10 -16 because of rounding errors in IEEE 754 double precision. The conclusion 
of Theorem 13.21 is employed to determine the time step at which we stop the calculation of the 
scheme 



5 Three-dimensional case 

We conclude the paper by a brief account of our numerical scheme in the three-dimensional case. 
Let = (0, Li) x (0,L 2 ) x (0,i 3 ) be a rectangular domain in R 3 , and put X = SI x S 2 . Let Mi, 
M 2 , M 3 , Mg, and be positive integers, and At be a positive number. We consider the grid 
points (tk,Xiji,£ mn ) given by 

A Xl = Li/Mi (i = 1, 2, 3), A9 = tt/M 9 , A<j) = 2tt/M , 
t k = kAt, Xiji = (iAxi, jAx 2 ,lAx 3 ), 9 m = mA9, 4> n = nA0, 
) = (sin 9 m cos <f> n , sm9 m sm<f) n , cos6> m ), 

where k, l,m,n € Z. Using I?, l m n as the value corresponding to I(tk,Xiji,^ mn ), we discretize 
the problem (fT]) for d = 3 as follows: 

A jk y jk+1 _i_ 7V- jk i fJ k 

& i,j,l,m,n & i,j,l,m,n ' A 1 i,j,l,m,n ' "i,j,l,m,m 

0<k< T/At - 1, (x ljh Un) G X, (24a) 
Io(xiji,£mn), (Xijl,£mn) G X, (24b) 

h(t kl x l3U Un), < k < T/At, (xiji^mn) e r_, (24c) 



i r +1 

-L i,j,l,m,r. 



Tk 

i,j,l,m,n 



c(Xiji) 



At 



i,j 7 l,m 7 n 



i,j,l,m,n 
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t = 50 [ns] 




0.0001 
0.0002 
0.0005 
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0.1 



1 = 100 [ns] 




0.0001 
0.0002 
0.0005 
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0.002 
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0.1 



(a) t = 50 (k = 500) 
t = 200 [ns] 

50 




0.0001 
0.0002 
0.0005 
0.001 
0.002 
0.005 
0.01 
0.1 



(b) t = 100 (k = 1000) 
t = 400 [ns] 




(c) t = 200 (k = 2000) (d) t = 400 (k = 4000) 

Figure 2: Contour of Integrated Light Intensity 
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Figure 3: Difference between Numerical Solutions and Steady State 



where q^^^ m ^ n = q(t kl x tJ i,£ m n), and the operators Aa, Sa, and are defined by 



i,j,L,m,n 

jk jk 

s- i-\-l,j,l,m,n i—l.j,l,m,n , i £ i 

Smn,l 2Ax |vron,l | 

jk jk 

j. i,j-\-l,l,m,n i.j — l,l.m,n i ^ 



e 



mn.3 " 



2Ax 2 

jk jk 

iij,l-\-l,m,n i,j,l — l,m,n 



rk 

i-\-l,j,l,m,n 


— 1T k 

i,j,l,m,n 


' i— l,j,/,m,n 








Tk 

i,j + l,/,m,n 


- 2I k , 

,m,n 


1 l,t,m,n 




2Ax 2 




Tk 

i,j,l-\-l,m,n 


-2I k ., 


4- 

1 i,j,l — l,m,n 


2Ax 3 



2Ax 3 

M e -1 

K Ali,j,i,m,n = H s {xni)A6A(i> P( X W> W.^sin^/^-^^. 

p=l t/=0 

Sufficient conditions for the stability and convergence are given as follows. 

Theorem 5.1. (Stability and positivity for d — 3) Suppose that the discretization parameters 
satisfy 

c+At c+At c+At 



Axi Ax2 Ax3 



< 1, 



ae 2 



■p(x;Z,?(B,<l>))BmB 



Acf> 2 



■ sup 



d 2 



, 6 , 



Ad 2 sup 
where we put 

£'(0,4>) — (sin 6 cos <f>, sin sine/), cos#) 
and the supremum is taken over (x,£) G X, < 9 < n, < <f> < 2tt in (I25b|) . Then we have 

Halloo < ||/o||oo + HAHoo + c + ||g||ooT, < k < T/ At. 
Moreover if q > 0, Jo > 0, and 7i > 0, we have i£ • j n > for all k, i, j, l,m,n. 



(25a) 
(25b) 
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Proof. We can show that the solution of satisfies 



{l + cAi(/i s +/i a )} 



jfe+i 



< \\r 



■ cAt j l m n \ 



-cAt\\q\ 



(26) 



by repeating the argument which lead to ([12")) in Theorem 12.21 Next we shall show that 



Mg-l M4-1 

A6A(j) ^ ^ pjXijllZmn, £,pu) Sin ^ < 1 + (^ijhCmn)GX, (27) 

which corresponds to (IT4l) . The proof proceeds as in Theorem 12.21 by replacing (fT2|) and (Tl4l) with 
(|2l)| and (dTJ), respectively. 

To complete the proof we derive (f2~T|) . Fix (a;^/, £ mn ) £ X and consider the functions 



:=p(^;;£mn,£'(M)), (28) 
G(6):= g{e,^)d<P-A<pY.9{0,^). (29) 

By the mean value theorem for integrals, there exists a point 6' £ (0, 7r) such that 

/ G(6>) sin dO = irG(6') sin 9'. (30) 

Applying Lemma [2~3l to the function /(0) = g(9', <fi), we find a point 0' € [0, 27r) such that 

G(d') = ^A^g^e 1 ^'). (31) 
We put dMI and into O to see that 

■A0 J] / g(e,<f> v )smede = —Acp 2 g H {e',cj)')sme'. 



JO 



<p) sin 1 



i/=0 



By the standard error estimate for the composite trapezoidal rule (see, e.g., [5] Theorem 12.1), we 
have 

^ g{0,4> v )sm6d6 = AO £ </>„) sin0 M - ^A0 2 — cf> u ) sini 



for some 6^ £ (0, 7r). Thus 



•7T />27T 

/ 

JO 



) sin( 



Me-l 

A0A0 ]T 5 (0 M ,0„) sin 0, 

!/=0 l_l = l 

Mi-1 



7T 

12 



c9 2 



^A0 2 ff ^(0',0')sin0'-^Afl 2 A^ £ ^^sinfl 



t/=0 



(32) 



ci 2 



Since 0) sin# is continuous on0<#<7r, O<0< 27r, the mean value theorem shows that 

there exists (0",</>") G [0,tt) x [0,2tt) such that 







A< ^ 51 q^9(0 A.) sine 



?=e, 



= 2tt -^9(6, 4>") shx9 



(33) 
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Noticing that 

pTT />27T p 

f f g(0, 4>) sin 6 d<f> d6 = / p(xiji; £ mn , £') dw^ = 1, 
Jo Jo Js 2 

we obtain 

Mg-l M^-l 

^ Yl 9(9», M sin 0^ < 1 + fi* 

fj_=l u =0 

from (0S1), and (l25bl) . □ 

Theorem 5.2. (Convergence for d = 3) Suppose besides the hypotheses of Theorem 15. II that 
(i) the mesh ratios At/Axi (i — 1,2,3) are kept fixed; (ii) the coefficients // a and /x s are bounded; 
and (iii) the solution I of the problem JlJ belongs to C 2 ([0, T) x (Jf ur_)J and is bounded together 
with its derivatives of the first and the second order. Then we have 

\\I(t k , ; -)-I k \\oo< C(At + A9 2 + A4> 2 ), 0<k< T/ At, 

where C is a positive number independent of At, AO, and A(j). 

The proof is completely similar to that of Theorem 12.41 and is omitted. 
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